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This review discusses detector physics and Monte Carlo techniques for cryogenic, radiation de- 
tectors that utihze combined phonon and ionization readout. A general review of cryogenic phonon 
and charge transport is provided along with specific details of the Cryogenic Dark Matter Search 
detector instrumentation. In particular this review covers quasidiffusive phonon transport, which 
includes phonon focusing, anharmonic decay and isotope scattering. The interaction of phonons in 
the detector surface is discussed along with the downconversion of phonons in superconducting films. 
The charge transport physics include a mass tensor which results from the crystal band structure 
and is modeled with a Herring Vogt transformation. Charge scattering processes involve the creation 
of Neganov-Luke phonons. Transition-edge-sensor (TES) simulations include a full electric circuit 
description and all thermal processes including Joule heating, cooling to the substrate and thermal 
difi'usion within the TES, the latter of which is necessary to model normal-superconducting phase 
separation. Relevant numerical constants are provided for these physical processes in germanium, 
silicon, aluminum and tungsten. Random number sampling methods including inverse cumulative 
distribution function (CDF) and rejection techniques are reviewed. To improve the efficiency of 
charge transport modeling, an additional second order inverse CDF method is developed here along 
with an efficient barycentric coordinate sampling method of electric fields. Results are provided in 
a manner that is convenient for use in Monte Carlo and references are provided for validation of 
these models. 
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I. INTRODUCTION 

Cryogenic radiation-detectors that utilize ionization, phonon and / or scintillation measurements 
are being used in a number of experiments. Both the Cryogenic Dark Matter Search (CDMS) [Il[2] 
and EDELWEISS 3] dark matter search utilize silicon and / or germanium targets to detect recoils 
of radiation in the target masses. A combination of ionization and phonon readout is used to 
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provide discrimination of gamma- and neutron-recoil types. The CRESST dark matter search 
utihzes CaW04 targets and readout scintiUation and phonon signal to discriminate between recoil 
types. The advantage of reading out both phonon and ionization (or scintillation) signals comes 
about from the differing ratios of ionization and phonon energy or scintillation and phonon energy 
created in electron- and nuclear-recoils in the detectors. The ratio of these two energies leads to a 
powerful discriminator for the experiment's desired recoil type. 

Both the ionization and phonon readout can be used to generate position estimators for the 
initial radiation interaction, leading to fiducial volume selection. In the ionization signal this is 
generally accomplished by instrumenting different parts of the detector with independent readout 
channels and vetoing events with large signal contribution outside of the desired fiducial volume. In 
the phonon signal it is generally required to measure the early, athermal component of the phonon 
signal which still retains a position dependent component. 

The physics required to accurately model these detectors is presented in this paper along with 
appropriate numerical tricks that are useful for an efficient detector Monte Carlo. This paper pro- 
ceeds with a review of radiation interactions, charge transport physics, phonon transport physics, 
instrumentation. Monte Carlo techniques and relevant physical constants are included where ap- 
propriate. 

This paper will focus on the use of silicon and germanium detector masses, both of which are group 
IV semiconductors. However there are other relevant materials in use such as calcium tungstate 
(CaW04) which leads to a small loss of generality. 

A. The CDMS Experiment and Detectors 

The Cryogenic Dark Matter Search pQ |2] utilizes silicon and germanium detectors to search for 
Weakly Interacting Massive Particle (WIMP) dark matter [U |5| candidates. The silicon or ger- 
manium nuclei provide a target mass for WIMP-nucleon interactions. Simultaneous measurement 
of both phonon energy and ionization energy provide a powerful discriminator between electron- 
recoil interactions and nuclear-recoil interactions. Background radiation primarily interacts through 
electron-recoils whereas a WIMP signal would interact through nuclear-recoils. The experiment is 
located in the Soudan Mine, MN, U.S.A. 

The most recent phase of the CDMS experiment has involved fabrication, testing and commission- 
ing of large, 3 inch diameter, 1 inch thick [100] germanium crystals. The CDMS-iZIP (interleaved 
Z-dependent Ionization and Phonon) detectors are 3 inches in diameter and 1 inch thick with a 
total mass of about 607 grams [5]. The iZIP detector utilizes both anode and cathode lines on 
the same side of the detector similar to a Micro-Strip Gas Chamber (MSGC) as shown in 

Figure [l] and [2j Unlike an MSGC however, there is a set of anode and cathode lines on both sides 
of the detector. This ionization channel design is used to veto events interacting near the detector 
surfaces. An amorphous silicon layer, deposited under the metal layers, increases the breakdown 
voltage of the detectors. The total iZIP aluminum coverage is ~4.8% active and ~1.5% passive per 
side. 



II. RADIATION MODELING 

When using a Monte Carlo of a detector, it is often helpful or necessary to have a numerical 
model of radiation interactions in the detector. Many readers will find it valuable to use separate 
modeling software such as GEANT4 [lU]. A brief description of these interactions follows. 



FIG. 1: (left) A CDMS "iZIP" detector with photolithographically defined phonon sensors. The crystal is 
3 inches in diameter and mounted in its copper housing. The top surface contains an outer, guard phonon 
sensor and three inner phonon sensors from which an event's position estimate can be made. The opposite 
face (not shown) has a similar channel design, but rotated 60 degrees. 

FIG. 2: (right) Close-up view of the iZIP phonon channel and ionization channel (thin lines in between 
the phonon sensors). The phonon channel is held at ground and the ionization channel is held at ~ ±2 V 
for the top (bottom) surfaces. 

Low energy gamma-rays (x-rays) predominantly interact via photoelectric absorption in which 
all of the gamma-ray energy is deposited in a single interaction location. High energy gamma-rays 
interact via Compton scattering in which some of the gamma-ray's initial energy is transferred to 
an electron and the gamma-ray continues along a different trajectory with reduced energy. The 
gamma-ray will generally continue to scatter until it leaves the detector volume or terminates with 
a photoelectric absorption. In silicon (germanium), for photon energies greater than 60 (160) keV, 
Compton scattering dominates [71 [TT] . 

Both of these electron interactions result in a high energy electron being produced which then un- 
dergoes a rapid cascade process resulting in a large number of electron-hole pairs [HIIT^] . This initial 
cascade process ceases around the scale of the mean electron-hole pair creation energy (Eeh.create) 
resulting in an expected number of electron-hole pair n^fi = / Ef,h,create- Due to correlations 
in the cascade process, the variance in the number of electron-hole pairs is reduced, relative to 
Poisson statistics, and given by tig^ = F x Ueh, where F is the Fano factor T3]. These high energy 
electron-hole pairs will then shed phonons until they reach the semiconductor gap energy -Egap 
which results in the fraction 1 — -EgapZ-Egh, create of energy in prompt phonons and the remainder in 
the electron-hole pair system. 

Neutrons that interact in the detector bulk will knock an ion out of its lattice site and displace it 
to some other location in the crystal. This high energy ion will interact with both the lattice ions 
and / or valence electrons, with competing cross sections, before reaching some other location in 
the crystal. Interactions with the lattice ions can be described to first approximation as Rutherford 
scattering |14' with differential energy loss per unit length —dE/dx '-^ where v is the ion's 
velocity ^15j. For interactions with valence electrons, the number of electron states which can be 
excited to an accessible state outside of the Fermi sphere scales like velocity v hence the differential 
energy loss per unit length scales like —dE/dx ~ v |15j . A description of this process is shown in 
Figure |3] There are important screening and velocity dependent cross sections in both energy loss 
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potential and hence can only change 
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Two hypothetical ion-electron interactions, shown 
before and after the scatter and in the crysal and 
ion reference frames. In the parallel interaction, the 
final electron velocity is too small to scatter to an 
unoccupied region outside of the Fermi sphere 
and is forbidden. In the anti-parallel process, the 
electron can be scattered to an unoccupied state 
outside of the Fermi sphere. 
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FIG. 3: Description of ion-electron potential interaction process which shows that for ion velocity much 
less than the Fermi velocity {v <^ Vf) the number of electron states that an ion can interact with scales like 

V. 



lengths P^SHTO] but to first approximation these equations show the velocity dependent competition 
in the two scattering rates. Whichever interaction occurs first, there is generally a cascade of 
many scattering events resulting in a larger amount of energy being deposited in the ion lattice 
compared to gamma-ray interactions. For historical reasons, the reduced amount of energy in 
the electron-hole (ionization) system, compared to an equal energy deposition by a gamma-ray, 
has brought about the term nuclear quenching to describe the reduced ionization signal. The 
details of the cascade and reduction in the number of electron-hole pairs is described by Lindhard 
theory [H] and given as / = kg{€)/{l + kg{e)) where e = ll.bEnZ^'' , k = O.lSSZ^/^A-i/^ ^nd 
g{e) — Ze^-^^ + 0.7e°'^ + e [50]. The recoil energy is En and given in units of kcV, Z is the atomic 
number and A is the atomic mass. 

Beta radiation represents another class of electron-recoil interactions. The attenuation lengths are 
much shorter, however, resulting in a class of events sometimes referred to as surface interactions. 
These surface interactions can result in signals that differ from bulk events of the same energy and 
electron- / nuclear-recoil type. For example, the events may be located in regions of large electric 
fringing fields, directing charges away from readout electrodes or near mounting hardware that 
absorbs scintillation light. 
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III. PHONON SIMULATION 



A. 



Introduction 



Phonons, in the context of this review, are quantized vibration modes that exist in periodic 
structures such as sihcon and germanium crystals. They are excitations which, in the lattice, and 
in materials sufhciently cooled that charge carriers are frozen out, mediate thermal transport. They 
are described by the Hamiltonian 



where m is the mass of the lattice atoms and u) is the frequency of oscillations about the center of 
the harmonic potential between an atom and its nearest neighbor atom ^21ri23j . 

In one dimension, the solution to the time-independent Schrodinger equation H \x) = E \x) yields 
the different oscillation modes at atom j 



where i = T here, the wave number fc„ = a is the lattice spacing and N are the number of 
lattice sites. 

In Monte Carlo, phonons are simply treated as non interacting particles with decay and mass 
defect scattering properties described in the remainder of this chapter. In general they have a 
nonlinear dispersion relationship; however, due to rapid down conversion to lower frequencies they 
spend most of their time in an energy region with linear dispersion relationship. Hence, a linear 
dispersion relation is sufhcient for phonon transport modeling in these detectors. 



Immediately after the recoil, a population of prompt phonons exist in the immediate vicinity of 
the interaction point. Particular details of the frequency distribution and mode population are not 
well known but we can make a few deductions, explained in more detail in the following sections, 
that will lead us an initial distribution for Monte Carlo. 

Anharmonic decay, due to nonlinear terms in the elastic coupling between adjacent lattice ions, 
causes the phonons to rapidly down convert into a lower frequency distribution. This process allows 
us to start a Monte Carlo with any high frequency distribution and details of the distribution 
will rapidly be lost; we use the Debye frequency as a naive starting point. Isotope scattering, 
which also occurs at a high rate for a high frequency phonons, causes the phonons to obtain 
their equilibrium mode density; we use the equilibrium mode density as a naive starting point. The 
approximations are good in the sense that the detector's phonon response is insensitive to variations 
in the distributions. 

They are valid since the detectors are large compared to the initial characteristic interaction 
lengths of phonons. Furthermore, later generations of phonons are ballistic and timing information 
in the measured phonon pulses is determined by the detector geometry and the loss rate of phonons 
at surfaces. 




(1) 



Xj 



(2) 



B. Prompt Phonon Distributions 
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C. Phase Velocities and Polarization Vectors 



The so called phase velocity surfaces represent the direction dependent phonon phase velocity 
Vp = Vp{S,(j)) . In general they are given by the eigenvalue Equation^ 



pw^e^ = X! X! ^t^-y^rKK Cr, (3) 

r \ (Tu / 

where p is the crystal's mass density, 
uj is the phonon frequency, 

is a component of the polarization vector e, 
c^iavT are components of the elastic constant tensor and, 
ka is a component of the phase velocity vector k [51]. 

Not all of the elastic constants are independent, reducing via a Voigt contraction [21] the number 
that we need to keep track off. Additionally, symmetries in a cubic crystal allow for further reduction 
in components and we can define three independent constants Cn = c^xxx = (^yyyy = Czzzzj 
C'i2 = Cxxyy = Cyyzz = c^zxx and C44 = Cxyxy = Cyzyz = Czxzx- This Contraction simphfies 
Equation [3^ significantly |22j and we are left solving matrix H for its eigenvectors and eigenvalues. 



^x^x C/^^ikyky -\- kzkz) (C'12 + C/^/^kxky 

{C12 + C44)kxky Cii{kxkx + kzkz) + Ciikyky ■■■ 

(C12 + Cii)kxkz (C12 + Cii)kykz 

(C12 + Ciijkxkz 

{Ci2+Cii)kykz 1 (4) 

C44{kxkx + kyky) + Ciikzkz 

The eigenvectors represent the three polarization vector directions and the three eigenvalues 
equal puj'^ for the longitudinal, slow-transverse and fast-transverse modes. The anisotropy in the 
cubic silicon crystal leads to the phase velocity surfaces being non-spherical (see Figure |4]). The 
phase velocities are used to determine both the group velocities (Section jlllD ) and also the isotope 
scattering rates. 



D. Group Velocities 



Phonon group velocities are found by solving 



dk 



(5) 



The slight lack of sphericity in the phase velocity surfaces (see Figure |4| has a very dramatic 
effect on the transverse phonon group velocities (25 -.28J (see Figure [s]). The longitudinal phonon's 
group velocity is only mildly affected. Energy is focused in the direction of heavy banding and 
leads to the term phonon focusing. The point density in the plots is misleading as the three modes 
are shown to be equally populated. Isotope scattering including anisotropic scattering rates (see 



Section III E ) leads to the phonon modes in silicon being populated as follows: slow-transverse 



(55%), fast-transverse (35%) and longitudinal (10%). 
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FIG. 4: Phase velocities in silicon. The distance from the origin represents the phase velocity speed (in 
m/s). The first three plots are views from the north pole. The fourth view containing all three surfaces is 
offset from the north pole. In the plots all modes are equally populated, which does not reflect the actual 
mode populations. 



E. Anisotropic Isotope Scattering 

Plionons scatter off mass defects in the crystal (see Figure |6]). Additionally, tliey can change 
modes. The bullc scattering rate is given by 

Ti = B,y%s% (6) 

where v is the phonon frequency and _B is a scattering rate constant E51 - I5T] (see Table |l| . The 
scattering rate for individual phonons is given by 

7-^^, (7) 

where vx' is the final state phonon frequency in Hz, e is the polarization vector, A represents the 
initial phonon and A' represents the outgoing phonon [28j [31] . It is the dot product in Equation [t] 
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FIG. 5: Group velocities in silicon. Energy is focused in the direction of heavy banding and leads to the 
term phonon focusing. The distance from the origin represents the speed that phonon energy is carried 
through the crystal (in m/s). In the plots all modes are equally populated, which does not reflect the actual 
mode populations. 



which allows mode mixing and the denominator which ensures the correct populations in the ratios. 
In silicon the populations when including anisotropic scattering rates are slow-transverse (55%), 
fast-transverse (35%) and longitudinal (10%). The standard treatment is to determine if a phonon 
isotope scatters via Tj and then determine its polarization and direction via 7. After the initial 
anharmonic decay has settled down, isotope scattering dominates. 

This process is unfortunately computationally expensive due to sample-rejection techniques [3 2) 
and after several iterations an isotropic scattering process can be used for individual phonons with 
little loss of modeling accuracy. 
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FIG. 6: Phonons isotope scatter off mass defects in the crystal. Equation [7] gives the individual plionon 
scatter rates. iTg is the group velocity and e is the polarization vector. 



F. Anharmonic Decay 



1. 



General Case 



Nonlinear terms in the elastic coupling constants cause a longitudinal phonon to down convert 
to two lower energy phonons (see Figure m). The bulk decay rate is given by 



where z/ is the phonon frequency and A is a decay rate constant PSI - IST] (see Table |T]). The decay rate 
for transverse phonons is negligible 33 . The three body problem requires that both energy and 
momentum are conserved. These conditions make an exact solution computationally prohibitive 
for large amounts of phonons. 



To allow computations to proceed in a finite time, an exact solution to anharmonic decay is 
abandoned and an isotropic approximation is used (the full anisotropic phase velocities and group 
velocities are still easily used for isotope scattering and phonon transport). The energy distribution 
calculations are still difficult but fortunately have already been carried out [311 133- Once the 
energies have been determined, calculating the resultant scattering angles based on energy and 
momentum conservation is fairly straightforward. Due to the different energy-momentum dispersion 
relations for the longitudinal and transverse phonons, there are two different decay branches, L — ^ 
L' + T' and L^T[+T!^. The L' +T' energy distribution is given by 



(8) 



2. Isotropic Approximation 



1 




(9) 
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FIG. 7: Longitudinal phonons decay due to nonlinear terms in the elastic coupling constants. tTg is the 
group velocity and e is the polarization vector. 



This approximation results in the outgoing phonons having an angular displacement from the 
initial phonon given by 



cos(6'l') 



2x 



(10) 



COS[t)T' 



l-a;2 + (52(i_a;)2 



2S{l-x) 



The L T[ + energy distribution is given 



^L^Ti+Ti - {A + B5x-Bx^f 



Cx{S - x) 



D 



6-. 



ix 



n 2 



(11) 



(12) 



where x — Sj^, 

si'"* SI 
2 2 ' 

A^l{l^S^)[l3 + \ + {l + d^){^ + ti)], 

B = P + \ + 2S^{-f + n), 

C = /3 + A + 2(7 + /x) and, 

D = (1-52)(2/? + 47 + A + 3m). 

The constants 7 and fi are the Lame constants and /3 and 7 are third order elastic constants in an 
isotropic model (there is additionally a third independent elastic constant a but it drops out of the 
equations). 
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This approximation results in the outgoing phonons having an angular displacement from the 



initial phonon given by Equations 13 and 14 



cos{erO = , (13) 

cos(^T,) = 2^(1 ^ (14) 

where the trivial substitution a; —> 1 — a: is made for the second phonon. With these closed-form 
energy densities and scattering angles, plots can be generated to aid understanding of these events 
(see Figures [S] and [9]) . 

These angles are relative the initial momentum vector k and need to be converted into the Monte 
Carlo coordinate system. Polar coordinates are useful and angles are provided in this system. In 
addition to rotation angles 0l' and 9t' (or Or^ and Ot^) an additional azimuth angle, relative to 

the initial momentum vector k, and in the isotropic approximation is randomly distributed from 
[0, 27r], is specified as 92tt- 

In terms of initial elevation and azimuth angles and 0, scattering angles Bj^^ and 0^^ and 
azimuth scattering angle ^27r, the final angles $1 and Oi that describe the phonon momentum 
vectors fci are 

f&i = arccos[— sin $ sin 6t^ cos 62-^ + cos $ cos 6t^] (15) 

and 

01 = G — arctan2[— sin 6x^ sin 6*2^ sin $, cos 9t^ — cos cos $1]. (16) 

The final angles $2 , ©2 for the other phonon momentum vector fc2 are found by replacing 6ti 
with and 027r with 27r — 02ir- 

G. Phonon Losses at Surfaces 

Eventually the phonons will interact at surfaces where they are instrumented, reflected back into 
the crystal, down converted to a lower energy or are lost to the environment. 

Phonons can reflect either specularly or diffusively from the surfaces. Specular reflection can occur 
on smooth, untreated semiconductor wafer surfaces. It is the simplest to describe with incident and 
reflected angles relative to the normal equal, 6i = Or- 

Diffusive scattering is also common on surfaces that have been damaged or roughened during 
fabrication. In the ideal case the scattering angle is described by Lambert's cosine law where 
the angular distribution scales like cos 6 where the angle is measured relative to the normal. 
Scattering surfaces that satisfy Lambert's cosine law scatter phonons isotropically regardless of 
their incident angle. Generally diffusive scattering has been found to be a good model for phonon- 
surface reflections, likely due to some small roughness in the surface. 

Phonons are strictly-speaking eigenstates of a Hamiltonian that describes an infinitely large, 
periodic lattice. This description necessarily breaks down at the detector surfaces and there is some 
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FIG. 8: Resultant energies in longitudinal phonon decay in germanium. The two branches L ^ L + T 
and L ^ T + T are shown; the plotted distribution is indicated in bold face in the legend. 




FIG. 9: Resultant angles in longitudinal phonon decay in germanium. The two branches L — >■ L + T 
and i — > T + T are shown; the plotted distribution is indicated in bold face in the legend. 
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probability that phonons down convert to lower energy daughters at the surface. The details of 
this process would be highly material / surface treatment dependent but the probability of this 
occurring for a particular phonon-surface interaction will be small for a high purity crystal. It is 
generally easiest to tune this probability by running numerous Monte Carlo to find the best value. 
In the CDMS detectors we have measured a loss of ~0.1% for each phonon-surface interaction 

It is the goal of the experimental setup to absorb the phonons into some sensor and provide instru- 
mentation into a data acquisition system. The details for the phonon-sensor interaction probability 
are again complicated and highly depend on the type of absorber and attachment / fabrication de- 
tails. Acoustic mismatch theory provides a good starting point for analytic calculations and can be 
performed over both normal and non-normal incidence angles. The relevance of these calculations 
can be lost however when detector to detector variations are considered. Additionally the angular 
dependence of such calculations can be washed out when integrating over a distribution of phonon 
incident angles and phonon energies. In practice it is usually again easiest to tune this probability 
by running numerous Monte Carlo and identifying the best fit value. In the CDMS detector there 
is additionally an amorphous silicon dielectric in between the crystal and thin metal films; we have 
empirically matched a phonon-aSi-aluminum interaction (details of the interaction are discussed in 
Section IV) probability of ^33% with the remaining 67% diffusively scattering back into the crys- 



tal [36ll37j. For any well designed detector, phonon absorption into the instrumentation sensors will 
dominate over other loss processes allowing the probability to be tuned by matching pulse decay 
times. 



H. Time Steps 



It would be grossly inefficient to run all phonons with the same time step. Therefore we gener- 
ate scattering times according to the distributions in Sections |III E| and |IIIF[ The scattering and 
decay probabilities go like P = 1 — exp(— </r) where r is a combination of isotope scattering and 
anharmonic decay rates as given by Matthiessen's rule 1/r = 1/Tiso + 1/Tan/i- 

This scattering time has to be compared with the time that the phonon will take to interact 
with a surface and the event that occurs soonest will be chosen for each individual phonon. If it 
is determined that the bulk interaction time is less than the surface interaction time, then it must 
be determined which event occurs based on their relative, frequency dependent rates. This can be 
done by drawing a uniform random number u and comparing the rate of the process in question to 
the total rate. For example, if u < {1 / Tank) / / Tank + 1/Tiso) then an anharmonic decay is selected. 



I. Random Number Sampling 



Only in rare circumstances will a uniform random number u be needed in Monte Carlo without 
some transformation. Often we are trying to draw the number x out of the probability distribution 
function (PDF) /. An efficient method for transforming u to the desired probability distribution 
function / is to integrate / to find the cumulative distribution function F = J f. The cumulative 
distribution function (CDF) F has the desirable property that it is bounded by [0, 1] as is u. The 
CDF F is then inverted and solved at u to determine f,x = F~^{u) [35] • 

As an example, we first considering the bulk interaction rate where the probability of having 
an interaction is P = 1 — exp(— t/r). After integration and inversion, the randomly generated 



15 



scattering time is given by 

^scatter = -Tln(M). (17) 

As a second example we can consider diffuse scattering ofl the detector walls. In a spherical 
coordinate system where 9 represents the angle from the surface normal, the PDF in this coordinate 
is / = sin(6'). This is also easily integrated and inverted to yield = arccos(u). Care must be 
taken in this example however since arccos is defined over the domain [-1,1] which modifies to 
F^^ = arccos(2M — 1). 

There are times when a PDF cannot be analytically integrated to yield a CDF or the CDF cannot 
be inverted. This is an unfortunate situation since an expensive rejection technique is required. 
This technique involves drawing a pair of uniform random numbers [ux,Uy] where min(.T) < < 
max(x) and min(/) < Uy < max(/). If Uy < fiu^) then is retained, otherwise is rejected 
and the process is repeated. The inefficiency of this method is related to the area coverage c = 
/ /(a;)da;/(max(/) x {ma.x{x) — min(a;))) and Ux will be successfully drawn in one of c attempts. 

The rejection method can be improved however for static distributions that do not change during 
the Monte Carlo run. An example includes diffusive scattering off of the side walls of the detector 
when using a spherical coordinate system. In this case, the sin(6') Jacobian results in the PDF 
/ = sm{9)^ for which F cannot be found analytically. The PDF can be integrated numerically to 
generate a CDF which is subsequently inverted. The process lacks a certain degree of elegance but 
is significantly more efficient than using a rejection method. 

J. Numerical Constants for Phonon Simulations 

Table |T] lists numerous constants that are used in the phonon simulations. They define the 
propagation dynamics, scattering and decay rates and energy carrier statistics. 

IV. QUASIPARTICLE DOWN CONVERSION 

Phonons have some probability of entering and interacting in the thin aluminum films that are 
patterned on the surface. These aluminum films make up both the phonon collecting films and 
ionization ground lines; it is interactions in the former that are measured in the phonon sensors. If 
the phonons have energy greater than or equal than twice the superconducting gap > 2ii^gap — A 
then Cooper pairs can be broken creating two quasiparticles. These quasiparticles will be of high 
kinetic energy i^kinctic = — A and contain some probability of scattering off phonons. These 
daughter phonons thereby introducing a population of down converted phonons back into the metal 
film. These phonons could break additional Cooper pairs in a cascade process that ceases when 
all phonons have energy below A or have a probability of being reintroduced back into the crystal. 
The key points in the cascade process are summarized in the following list [551 - HT| . 

1. Quasiparticle recombination lifetimes are long compared to quasiparticle decay and quasipar- 
ticle absorption into the aluminum films and therefore the recombination processes can be 
ignored. 

2. Quasiparticle decay via absorption of a phonon is suppressed at low temperatures due to a 
phonon density of states term n{fl), where Q is the phonon energy, in the Green's function 
and therefore can be ignored. 
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TABLE I: Numerical Constants for Phonon Simulations 



Parameter Descrip- 
tion 


Symbol 


Units 


Silicon 


Germanium 


Reference 


Isotope Scatter Rate 


(B) 




2.43 X 10"*^ 


3.67 X 10"*^ 


m 


Anharmonic Decay 
Rate 


(A) 


[A 


7.41 X 10"^'' 


6.43 X 10"^^ 


m 


Longitudinal 
Velocity 




[m/s] 


9000 


5310 


m 


Transverse Velocity 


M 


[m/s] 


5400 


3250 


m 


Decay Rate Constant 


(P) 




-0.429 


-0.732 


m 


Decay Rate Constant 


(7) 




-0.945 


-0.708 


m 


Decay Rate Constant 


(A) 




0.524 


0.376 


m 


Decay Rate Constant 


(^) 




0.680 


0.561 


m 


Density 


(P) 


[g/cm=*] 


2.33 


5.32 


m 


Decaying Ratio 

L—¥L-\-T 
All Lortqitudinal Decays 






0.204 


0.260 


m 


Lattice Constant 


(Cii) 




1.66 


1.29 


[22] 


Lattice Constant 


(C12) 




0.64 


0.48 


[22] 


Lattice Constant 


(C44) 




0.80 


0.67 


[22] 


Debye Frequency 




[1/s] 


1.5 X 10" 


0.864 X 10" 


[22] 


Electron-Hole 
Energy 




[eV] 


1.17 


0.75 


m 


Mean Electron-Hole 
Creation Energy 




[eV] 


3.81 


2.96 





3. Quasiparticle decay via emission of a phonon results in phonons with an energy distribution 
given by P^{fl) = fl'^p{E — ft) (^l — ^(^^-n) ^' where E is the quasiparticle energy and the 
quasiparticle density of states at temperature T—0 goes like p{E) — and < < i?. 

4. Phonons break Cooper pairs producing quasiparticles with an energy distribution given by 
Pqi;p{E) = (1 + Ei^-E) ) Pi^)pi^ - ■^)' where < E <Vl - E. The phonon is completely 
absorbed so that the second quasiparticle has energy E2 = ^ — E . 

5. Phonons are lost to the crystal if they reach the aluminum / crystal interface. 

A. Monte Carlo process ordering 

In the physical processes can be ordered as follows 

1. If the phonon energy is sufficient f2 > 2 x £^gap, a quasiparticle pair is created with the 
distribution Pqp(E) previously described. 
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2. If there is a quasiparticle with energy i? > 3 x -Egap then the quasiparticle emits a phonon with 
energy distribution P^(fl). Quasiparticles with energy i? < 3 x -Egap would shed a phonon 
with E < 2 X -Egap and therefore provide an endpoint for quasiparticle generation. They may 
be removed from the Monte Carlo. 

3. The probability of a phonon escaping the crystal is a function of the distance to reach the alu- 
minum / crystal interface and the phonon / Cooper pair interaction length. Given the large 
number of phonons it is generally not necessary to track this process in detail and instead a 
simple model is sufficient. On average, phonons are assumed to populate the center of the 
aluminum film (z — Zai/2), where is the aluminum thickness, and the phonons have 1/2 
probability of traveling upwards and 1/2 probability of traveling downwards. For the down- 
ward going phonons there is an exp{-~{2 x 1^i/2)/Lq) probability of reaching the aluminum / 
crystal interface before scattering, where Lq ~ 720 nm is a characteristic phonon interaction 
length [3T]. The factor of 2x is provided to integrate over different phonon incidence angles. 
The factor /ai/2 is replaced by Zai x 3/2 for upward going phonons. 

4. If the phonon has not been removed from Monte Carlo in step 2 or reintroduced into the 
crystal in step 3 then the process repeats at step 1. 

V. CHARGE MONTE CARLO 
A. Introduction 

Accurate modeling of charge propagation is included in Monte Carlo for numerous reasons. First, 
the ionization signal, compared to the phonon signal, provides a discriminator between electron- 
recoil and nuclear-recoil events in the silicon and germanium detectors. Second, electron transport 
is described by a mass tensor, leading to electron transport which contains components oblique 
to the applied field. This description is necessary to explain and interpret signals in the primary 
and guard- ring ionization channels, which function as a fiducial volume cut |42( I43j . Third, for 
electron recoils in the germanium bulk, charges drifting through the detector produce a population 
of Luke phonons which contribute 56% of the total phonon signal at ±3 V bias. This fraction 
is understood by considering that for every -Ech. create of gamma energy, an electron-hole pair is 
created which contributes eVbias of phonon energy; eVbias/(E'ch, create — E^gap + eH)ias) — 56% is the 
contribution of Luke phonons to the total phonon signal. These phonons' spatial, time, energy and 
emitted-direction distributions should therefore be properly modeled in Monte Carlo of the detector 
response. Fourth, phonons created during electron-hole recombination at the surfaces contribute 
13% of the total phonon signal but in a low frequency, ballistic regime that is used to provide a 
surface-event discriminator. This fraction is understood by considering that for every E^eh, create of 
gamma energy, Egap of phonon energy is released at the surface; eVgap/(Ech,croate + eVbias) = 13% 
is the contribution of electron-hole gap energy phonons to the total phonon signal. These phonons 
also need to be properly modeled in a Monte Carlo. 

Germanium has an anisotropic band structure described schematically in Figure [TO] and shows 
energy band structure in which the hole ground state is situated in the F band's [000] direction 
and the electron ground state is in the L-band [111] direction. Hole propagation dynamics are 
relatively simple due to propagation in the F band and the isotropic energy-momentum dispersion 

relationship e( fc ) = h^k^ /2m. Electron propagation dynamics are significantly more complicated 
due to the band structure and anisotropic energy-momentum relationship. At low fields and low 
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FIG. 10: Energy band structure of germanium, showing the L-valleys at (111), the F valley at (000), and 
the X- valley at (100). Symmetry results in 8 total L-valleys and 6 total X- valleys. 



temperatures, electrons are unable to reach sufficient energy to propagate in the F or X-bands, and 
are not considered necessary to consider in Monte Carlo. The electron energy-momentum dispersion 
relationship is given hy e{k) = x (fc|/m|| + fc^/mj^), where the longitudinal and transverse 

mass ratio m\\/m± ~ 19.5. 

This chapter will proceed by describing hole propagation and scattering, electron propagation and 
scattering utilizing a Herring- Vogt transformation and finally electron-hole recombination. Higher 
order mass terms and scattering processes which occur at high electric fields are discussed elsewhere 
in the literature 1441. 



B. Holes: propagation and scattering with isotropic bands and isotropic phonon velocity 



Hole propagation dynamics are described by momentum evolution in an electric field 



and propagation in position space it = hk /rric, where rric is the effective carrier mass. 

Charge carriers cannot accelerate indefinitely however, and the shedding of Neganov-Luke 
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FIG. 11: Charge carrier with initial wave vector k and final wavevector k' scattering off of lattice at angle 
(p with respect to k and emitting a phonon with wavevector q at angle with respect to k, where the vector 
momenta sum as shown on right. 



phonons |151 |3B] limits their speed to around the longitudinal phonon phase velocity v^. As de- 
scribed in Figure [TT] chargCj^honon scattering is an elastic processes, conserving both momentum 

k — k' = ~<t (where k and fc' are the initial and final hole momentum vectors and ~^ is the phonon 
momentum vector) and energy e — e' = hu (where e and e' are the initial and final hole energies 
and huj is the phonon energy). The phonon energy- momentum dispersion relationship is given by 

^ = v^q. Due to the low carrier energy, Umklapp processes |21) in which k — k' = l} + where 
G is a reciprocal lattice vector, are suppressed. 

Energy-momentum conservation coupled with the previous dispersion relationships leads to the 
final states fc'^ = k^ + — 2kqcos{9) and q — 2{kcos9 — ki,) and 



fc2-2fc,(fccos0-A:,)-2(fccos(0)-fc,)2 

cos(0) = , , (18) 

fc^(fc2-4fcs(fccos(6')-fc«) 

where 9 is the angular displacement between k and if , (p is the angular displacement between k 

and fc', and fc^ is defined as fc^ = mvL/h. If we can determine a scattering rate r and phonon 
angular displacement 0, then we can use these formulae to find the final states. 

Fermi's Golden Rule provides the transition probability per unit time per unit energy as 



Pk. 



2tiV 



k'±q 



< ± ~t\H\Y >' 6{E- E' T 



(19) 



For phonon emission processes. 



±~^\H\k' > 



q{iiq + 1) 



(20) 



2VpvL 

where C is the deformation potential constant and Uq is the phonon occupation number given by 



'"^q — ^1 (e!^^^^'^ ~ l)- ^ characteristic length can be defined as 
probability can be integrated over 9 and E' to obtain a scattering rate 



fe4 



The transition 
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l/r = 



3 ^0 fci 

The angular distribution then follows to be 



1 - 



(21) 



P(k,e)d0 =''^(Jl\ ( cos(9) - ^\ smOdO (22) 
h \kLj \ k J 

where < < arccos(fci/fc) < 7r/2. 

The phonon scatter azimuthal rotation angle is uniformly distributed about < < 27r. The 
charge carrier azimuthal rotation angle is required io he — -k -\- "d . 

It is critical that the time steps in the Monte Carlo are sufficiently small that the scattering 
rates are approximately constant during a time step. A method for ensuring this is discussed in 
Sections EH and IVEl 



C. Electrons: propagation and scattering with anisotropic bands and isotropic phonon 

velocity 

Electron propagation and scattering is complicated by the anisotropy of the electron bands 
but can be simplified by performing first a transformation into a space defined by the vectors 
(ej|, e_Li, ej_2) (where ejl is aligned with [111] and the other two are perpendicular) and then a 
Herring- Vogt transformation into a space where the electron bands are isotropic [12 HZj- The 
Herring- Vogt transformation is non-unitary and in the (ejl, eXt, ej^) space is given by 



V "'II 





\ 





/ 'Tflc 





V 





/ nic 



but the speed of sound vl remains unchanged and isotropic. 

In this space, k* — Tk, e — ^^^^ , — eE , and the effective mass is given by S/rric = 

l/my + 2/m±. The change in velocity, in position space, is found by a back transform, incorporating 
both the mass and momentum transforms, and is given by u = —Tk*. 

After the Herring- Vogt transform is used to find the electric-field induced acceleration, the elec- 
trons shed phonons via the same prescription given to the holes. The phonon and electron momen- 
tum is found first in the Herring- Vogt space and the back transform is applied to return to 
position space. The only additional concern is correctly back transforming the phonon momentum 
due to the non-unitarity of the Herring- Vogt transform. To handle this, we maintain the phonon 
momentum magnitude (ie, conserve energy), but use the back transform to find the correct 
angular distribution. 
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D. Charge Time Steps, First Order 



Determining an efficient time step size, dt for charge transport is complicated since the scattering 
time T varies as the charge carrier accelerates. In this charge transport model it was decided to use 
sufficiently small dt such that r is relatively constant at each iteration. The requirement, sufficiently 
small, in practice means that energy is conserved in the transport process and the following is a 
detailed description of how this requirement is implemented. 

By observing the charge Monte Carlo over numerous field strengths and in germanium, it was 
observed that scattering limits the maximum possible charge momentum magnitude to about 



'^max,el 



13 X kr 



1/3 



(24) 



kmax,h = X kL\'^\^^^ (25) 

for electrons and holes respectively. These momenta k are then used to determine stepping times 

which conserve energy; the shed Luke phonon energy must equal the change in potential energy 
eAV at the carrier drifts. Again running numerous Monte Carlo it is determined that a stepping 
time of dt = t/2 is sufficiently small to conserve energy. Larger stepping times will result in a 
deficit of Neganov-Luke phonons being created. 



E. Charge Time Steps, Second Order 



The earlier described first order method can be improved upon by developing a second order 
method. The challenge is to efficiently and accurately determine the time until a charge sheds a 
Neganov-Luke phonon, which is challenging due to the changing interaction time, as the charge is 
accelerated by the field. An inverse CDF technique would be advantageous and one is developed 
here that adapts to the changing interaction time. 

This is done by sampling the scattering rate at two diflferent times tq and ti . We start with the 
differential equation ^ = —aoN and expand to next order in time 

dN , 

— = (-ao - ait)N. (26) 

Integrating, we can obtain 

InN = -{aot + ait^l2). (27) 

This continues with the standard technique of solving for the CDF and inverting to obtain a^t^ /2 + 

Got + ln(u) = which can be solved for the scattering time t = 2aiirm ^ positive root 

is retained as physical which provides a scattering time of 
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This is completed by recognizing that 



ao 



dN 
~dt 



(29) 



and 



ai 



1 / dN 




dN 




ti — to \ dt 




~dt 





to 



(rr^-To-i). (30) 



There is a maximum samphng time step {ti — tp) that can be used before the Hnear interpolation 
of scattering rates is inaccurate. As in the first order case, this results in lack of energy conservation. 
It is found that the electron sampling time step can be a factor of ^15 greater than the time step 
shown in Equation [24] and the hole sampling time step by a factor of ~20 greater than the time step 
in Equation |25[ Given these sampling time steps, most Neganov-Luke phonons are produced in a 
time t < ti, hence this method is much more efficient than the first order method. The efficiency 
of this second order method also implies that pursuit of higher order methods will not yield much 
additional improvement in computational efficiency. 

This method also couples well to a second order spatial transport method. The velocity form of 
the Verlct algorithm [48, 49^ is convenient and given by a description that should look familiar. 

1. x{t + At) = x{t) + v{t)At + \a{t)At^ 

2. v(t + At) = v{t) + i(a(t) + a{t + At))At 

This can be easily modified to incorporate the second order inverse CDF sampling method via 
the following procedure: 

1. Make a guess for the step size At, which we will call Atg 

2. x{t + Ata) = x{t) + v{t)Ato + ^ait)Atl 

3. To = ro(u(t)) 

4. v{t + Ato) = v{t) + i(a(i) + a{t + At„))Ato 

5. Tl = Tl{v{t + Ato)) 

6. From second order inverse CDF method, determine the randomly distributed scattering time 
Ati 

7. If Ati < Ato 

(a) x{t + Ati) = x[t) + v{t)Ati + \a{t)Atl 

(b) T2 ^ Ti{v{t + Ati)) 

(c) v{t + Ati) = v{t) + i(a(t) + a{t + Ati))Ati 

(d) Save T2 and a{t + Ati) for use in next iteration 



8. 



Else, save ti and a[t + Ato) for use in next iteration 
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Since holes are described by a scalar mass, this procedure is straightforward. There is however 
a slight modification to this procedure that is useful for electrons. Due to the use of the Herring 
Vogt transform, it is generally easier to keep track of momentum in the Herring Vogt space, fc*, 
rather than velocity. We also identify that the acceleration is given by a = X^^^^XE, where the 
transform X from the cartesian space to the space defined by basis vectors (e|, e±i, ej^) is shown 
explicitly. 

F. Select constants for charge Monte Carlo 



TABLE II: Physical constants for Si and Ge crystals. The isotropic hole effective mass m^, and the 
anisotropic electron effective masses my and m± are || and _L, respectively, to the conduction valley axes, 
and conductivity effective mass 3/mc = i/m\i + 2/m_i_. The incident energy per final electron-hole pair is 
eeh, vl the speed of sound, and lo = 7rfi'*p/(2m^H^) is the characteristic range for carrier scattering where 
Si (from [13]) or H At (fit to data (50]) is the deformation potential. 





Silicon 


Germanium 




Electrons 


Holes 


Electrons 


Holes 


mh/rue 




0.5 




0.35 


m|| /rUe 


0.91 




1.58 




mj_ / me 


0.19 




0.081 




mc/me 


0.26 




0.12 




Vl (km/s) 


9.0 


5.4 


P (g/cm^) 


2.335 


5.323 


Hi (eV) 


9.0 


5.0 


11.0 


4.6 


Hfit (eV) 






11.0 


3.4 


lo (/xm) 


16.9 


7.5 


257 


108 



VI. ELECTRIC-FIELD LOOKUP 

A. Electric-Field Lookup from Triangulated Mesh 

A numerical electric field model is necessary for the charge transport described in Section |Vj 
The simplest model is a constant, longitudinally directed field. However it may be desirable to 
include fringing fields and details from the electrode structure. A more accurate model will utilize 
a triangulated mesh. A 3-d mesh contains nodes, with each mesh node containing an associated 
electric potential V. At points other than a mesh node, the potential V must be interpolated. The 
MATLAB programming language offers a few options for this interpolation, the fastest of which, 
utilizing a barycentric-coordinates linear interpolation via the TriScatteredlnterp class. The Com- 
putational Geometry Algorithms Library (CGAL) is available for C++ [51] though this paper will be 
presented for a MATLAB implementation. Furthermore, the barycentric transformation involves 
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a linear transformation which imphes that the electric field is constant within a triangulation, a 
property which can be exploited to speed up computation. 

MATLAB's method of looking up the potential, while very convenient, is not efficient considering 
the number of repeated field queries that occur before the carrier has moved to a location with a 
differing field. Efficiency can be improved by exploiting the fact that a charge remains within 
its triangle for numerous iterations and that the field is constant with the triangulation. On the 
contrary, MATLAB solves for the potential at every iteration in its lookup procedure, which is a 
bit slow. These repeated searches can be avoided but at the expense of significant code complexity. 
The effort is justified however as charge transport imposes a dominant computational expense in 
Monte Carlo. 



B. Barycentric Coordinates 



Given a triangulation (the mesh is made of tetrahedra in 3-space but the term triangulation 
often persists) with four node points ri,r2,r3, and r4, the arbitrary point r can be described by 
the barycentric coordinates Ai,A2,A3, and A4 where 

r = Airi + A2r2 + A3r3 + A4r4. (31) 
The additional constraint is imposed that 



Ai + A2 + A3 + A4 = 1. (32) 
The barycentric coordinates A become more intuitive when thought of as area (volume in 3-d) 



coordinates (see Figure 12). In this paradigm, consider the 2-d node points ri,r2, and along 
with the probe point r. To start with, let's normalize the area enclosed by ri, r2, and to 1 (this 
normalization is equivalent to the constraint Ai = 1). Then we can consider the three different 
areas enclosed by 1) r2,r3, and r (ai), 2) ri,r3, and r (02) and 3) ri,r2, and r {a^). It turns out 
that these areas (ai, 02 and a^) are identically equal to the barycentric coordinates Ai,A2, and 
A3, providing a quick and intuitive interpretation of the barycentric coordinates. The process and 
interpretation is the same in 3-d when volume is substituted for area. It is not actually recommended 
to calculate A through this procedure but to instead follow the procedure in Section |VI C[ 



C. Barycentric Coordinate Formulae 



In this section we derive formulae useful for solving the barycentric coordinates and electric- 
potential. We start again with the definitions given by equations [31] and 32 After separating 



equation 31 into the x, y, and z components we can solve for the A through the following linear 
procedures and the formula is written out explicitly below. 



T : 




X4 X3 — X4 

2/4 2/3 - 2/4 

Z4 Z-i - Z4 



(33) 



FIG. 12: Mesh node points ri,r2, and ra along with the probe point r. The areas ai,02 and 03 are 
identically equal to the barycentric coordinates Ai, A2 and A3. 




X — X4 

T-\r-n)=T-' I | , 

Z — Zi 



where A4 = 1 — (Ai + A2 + A3). Explicitly we can write out 



(34) 



det(r) 



{Z2 - Zi){xz - Xi) - (Z3 - Zi){x2 - .T4) 

{x2 - X4)(2/3 - 2/4) - (a;3 - 2:4) (2/2 - 2/4) 
(2/3 - yi){zi - Z4) - {yi " 2/4) (23 - Z4) 

{Z3 - Z4){Xi - X4) - {Zi - Z4){X3 - X4) 

{xs - X4){yi - 2/4) - - 0:4) (2/3 - 2/4) 
(2/1 - 2/4) (^2 - Z4) - (2/2 - 2/4) (^1 - Z4) 

{Zi - Z4){X2 - X4) - {Z2 - Z4){X1 - X4) 

\ - X4)(2/2 - 2/4) - (2:2 - a;4)(2/i - 2/4) y 



where ciet(T) = ... 

{Xi - X4) X [(2/2 - y4){Z3 - Z4) - (2/3 - yi){Z2 - Zi)] + ... 
{X2 - X4) X [(2/3 - y4){zi - Z4) - (2/1 - 2/4)(z3 - Z4)] + ... 
(2:3 - X4) X [(2/1 - y4)(22 - Z4) - (2/2 - 2/2) (^1 - Z4)]. 



(35) 
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The potential is simple to solve for in the barycentric coordinate system and equal to V = 
XiVi + X2V2 + A3V3 + X4V4. This procedure is more obvious when one considers the connection to 
area coordinates. 



D. Barycentric Coordinate Procedures and Shortcuts 

The procedure for finding which tctrahcdra a charge resides is not unique and a canned MATLAB 
procedure (or CGAL library in C++) can be utilized for this step. After the tetrahedra in which the 
carrier resides is determined, the electric field, E = ~W is computed. This can be performed by 
probing the potential at four locations (r, r + Sx, r + Sy and r + Sz) and computing gradients. The 
drawback is that this procedure requires conversion to barycentric coordinates and electric potential 
lookup for three additional points. These steps can be eliminated with some simple derivations that 
are outlined below. Most of these steps provide a conceptual framework and only the last step is 
actually computed. 

First we consider two points r' and r" = r' + 6x and find their associate barycentric coordinates 
A' and A". 



A'. 






— Xi 


A'2 




y' 










- Zi 



(36) 



and 



/ ^1 ^ 




/ 


x" 


-Xi\ 


A^' 






y" 








\ 


z" 


- Z4 J 



(37) 



where, as always, A4 is constrained by ^A = 1. Next we solve for the potentials V and V" 



V' = X[Vi + A^V^2 + X'sVs + X'^Vi = X'^ 



(38) 



and 



V" = X'lVi + X'^V2 + X'^Vs + X'^Vi = A"V, 



(39) 



where V= (Vi, 1^2, 1^3, V^4). 

It follows that = {V — V")/5x but another route is preferable. We can define ^Ai^2,3 = 

f Sx\ 

A'1,2,3 — A"i,2,3 = —T~^ 

V / 

The fourth element in SX is determined by remembering that ^ A' = ^ A" = 1 it follows that 
^ (5A = and therefore SX4 = —{SXi + SX2 + SXs). We can define a matrix related to but of 
size 4 X 3 to make future calculations easier, 
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T ^2,3 




T ^3,2 











(40) 



1,3 



We have now obtained a convenient and efficient method for calculating Ex 



Sx 



W6X 



Sx 





1 Si] 






■ 

















(41) 



(42) 



and similarly for Ey and Ez ■ Equations 40 and 42 are used in Monte Carlo and the fields are reused 
until the charge leaves the triangle bounSaries. 



E. Determining if a Charge Leaves a Triangle 



The assumption that is exploited in these computations is that the charge remains in a particular 
triangle for many iterations. This is true, especially in the detector bulk where mesh nodes are 
sparse, but eventually the charge will enter a different triangle. With the concept of area com- 
ponents, it is clear that one of the barycentric coordinates A will be zero when the charge is at 
a surface. By continuation, the same A will become negative when the charge leaves the triangle. 
This can be tested efficiently via the linear calculations already derived in Equation [34] When a 
charge enters a new triangle then T, T and E must be recalculated and stored to memory. 



VII. TRANSITION EDGE SENSOR SIMULATIONS 



A. Introduction 



Energy collected in the aluminum fins diffuses through the aluminum fins and is collected into 
the tungsten Transition Edge Sensors (TESs) [S51 with an efficiency of ~10%. The aluminum 
quasiparticle fins combined with the tungsten TESs together are referred to as Quasiparticle-trap- 
assisted Electrothermal-feedback Transition-edge-sensors (QET) [S3]. The quasiparticle downcon- 
version process was described previously in Section |IV| and here we focus on modeling of the TES 
sensor, which provides the electrical signal which is readout in the experiment. Compared to the 
phonon simulations, the TES simulations are a bit simpler; most of the electrical and thermal 
processes that make up the simulation are likely to be more obvious to the reader. 

There are, however, tricky numerical issues in solving for the TES voltages. A description of 
the process follows, and a flowchart (see Figure 13) helps in understanding it. First, the biasing 
conditions have to be determined. When a TES is run as a phonon sensor (they are also found in 
many x-ray detectors), the bias voltage is held constant. However, current -voltage (IV) character- 
ization curves are often studied in which the bias current is swept through a range of values. The 
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biasing circuitry and TES resistance R = R{T,I), where T is temperature and / is current across 
some physical extent, are modeled. The TES electron temperature is affected by three processes: 
1) Joule heating warms the TES; 2) Conduction into the substrate cools it; 3) Diffusion through 
the TES spreads heat. It is of course the goal to run the simulation both quickly and accurately. 



B. Electrical Circuit Modeling 



The TESs are relatively simple to model as each can be thought of as a one dimensional object 
(see Figure 16). Generally, a TES is modeled with a minimum of two nodes and as many as ^100 
if it is desired to reproduce normal-superconducting phase separation f55H57j . The functional form 
of the resistance at each node has generally been given by 



Ri,j — 



Rn 



— Rn 



1 + tanh 



(43) 



for li j small and Ri j 



Rrnax, toT li^j large; where Rmax is the TES's resistance in the normal state. 



Rmin is the superconducting resistance (numerically it is best to not set Rmin to zero but rather 
some small value ~ 10~* Tc is the midpoint temperature of the superconducting transition, T^, 
is the temperature width of the superconducting transition, Ic is the critical current (see below), 
and Use ~ 2/3 is motivated by [551 [55] . 

The critical current is provided by Ginzburg-Landau theory and near Tc is equal to Ic = 

3/2 



3.52 



where /cb is the Boltzmann constant, and Cn is the heat capacity. 



Special care needs to be taken when calculating the critical current when modeling TESs that are 
wired in parallel. Implicit in the Ginzburg-Landau equation is that one thin film superconductor 
contains the current while in CDMS detectors, the TES is distributed over tites TESs. This cause 
the heat capacity C„ to scale down by l/nxES and Ic scales down by l/.y?ixES for each of the 
TESs. There is an offsetting factor when the currents from all the nxES TES which causes Ic to 
scale up by tites- The end result is that Ic for the entire collection of tites parallel TESs scales 
like i/riTES- 

The requirement that the resistance R = Rmax at large current / is required for IV curves that 
sweep from the superconducting to normal states. There is no precise definition of large /, but it 
is chosen to be much higher than the quiescent operating current while low enough that it is not 
too much higher than the normal to superconducting transition current. A plot of i? = R{T,I), 



dR/dT X T/R and P = dR/dl x I/R is shown in Figure 14 



Modeling of the TES and biasing circuitry (see Figure 15 ) is complicated by our desire to be able 
to create IV curves in which the bias voltage iVbias) is swept from zero to some large value and 
back to zero again. This difficulty comes about from the current dependence on resistance and also 
due to the sharp change in current that occurs during the transition from the superconducting to 
normal states. This difficulty can be overcome if we use a damping factor to prevent rapid resistance 
changes. To determine the TES resistance at step n we let 



Rn — 



iC-l)Rn-l+R{T,I) 

c 



(44) 



where is a stabilization constant whose value ~ 5 is found by trial and error. Picking a large value 
makes current transitions numerically stable, but requires slower step sizes dt to retain accuracy. 
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<: 



Determine bias 



voltage 

■ <: 



Find TES voltage 
assuming Vbias = 1 










Find current through 
TES and all other 
circuit elements 






Correct TES voltage 
bias and find voltage 
and currents 



I = ( U + U-1 ) /"2| — I 



Determine net power 
and temperature 



Rn=((g-1)Rn-1 + R(T,l))/g[ — I 





FIG. 13: TES simulation flowchart. 
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FIG. 14: Surface and contour plots of i? = R{T,I), a = ^| = a{T,I) and /3 = = I3{R,I) for 

a high-Tc, inner iZIP channel. The colors in the surface plot indicate the value of resistance, alpha and 
beta with blue representing and red the highest value in the figure. The contour plots show the same 
information but over a limited current and temperature region. The black dot indicates a nominal bias 
region, which will affect noise and pulse shape after a radiation interaction in the detector. The gradient 
in resistance and temperature is generally along the temperature direction, whereas for /3 it is in a mixed 
— T + / direction. 



Next, we solve for the voltages at each TES node (to be discussed below) and determine the 
resulting voltages and currents. We also model the inductor as discussed below. It is likely that 
the currents that we solve for differ from those that we originally used to compute R — R{T, I). I 
keep evaluating RmiT, Im), V and / letting 

Im = ^ (45) 

until 



(46) 
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where r; ~ 10 ^ is satisfied for all nodes, ensuring a self consistent answer. I have enumerated the 
individual steps m here to prevent confusion with the steps n that involve simulation time steps dt 



(again, see Figure 13 1 



The rest of the circuit, including the inductor, is modeled by finding the macroscopic current 
through the TES /„. A simple circuit analysis reveals that it is given by 



" RbiasRshunt -L — ' 

D, 4_ R , ' ■'■ ^true ' -I- ''parasitic ' 



1. TES Voltage Modeling 

In the simulation, temperature, and therefore resistance, is modeled at each node. However for 
the purpose of modeling the electrical circuit the interconnects need to have a resistance which is 
given by averaging the resistance of the two nodes that they connect. 

We will now set up the system of equations GV = B where the conductance matrix, G, is a 
sparse, square matrix and i? is a vector that describes the electrical boundary conditions . The 
form of the conductance matrix G is simplified significantly if we treat the individual TESs as one 
dimensional object. In general (regardless of the 1-d simplification) nodes either are voltage biased 
(which includes grounding) or are not. Those that are not biased (or grounded) impose current 
conservation requirements. From Kirchhoff's current law, the equation to be satisfied for nodes 
that have four adjacent nodes is 



^Icft ~ ^ _|_ bright — Vtop — Vbottom — V 

R\clt bright -Rtop ^bottom 



= 0. (48) 



Nodes without four adjacent nodes will result in the appropriate terms being removed, for example 
a node with a top, left and right neighbors is given by 

Meft^^KigM^^^W^^O. (49) 
^Icft bright -Rtop 

The result on the G matrix is straightforward. The diagonal elements are given by 

G.. = - (50) 

^Icft bright 

The non-diagonal elements that couple adjacent nodes are given by 

Pterins which couple adjacent nodes — • 

(51) 

-n^adjacent 

Nodes that have their voltage fixed either through voltage biasing or grounding satisfy the equa- 
tion V = Vbias- They are given the coupling terms 

G,, = 1, (52) 
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the other terms in the row are 

G,, = 0, (53) 

where i ^ j. 

Now we can set up B and complete Kirchhoff's current law. For the rows in G that describe 
current conservation, 

B, = 0. (54) 
Whereas the rows in G that describe voltage biasing, 

B^ = l/bias. (55) 

We can then solve 

V^G-^B, (56) 

to find the voltages at all TES nodes. 

C. Thermal Processes 

There are many thermal processes which affect the TES temperature: the addition of heat from 
phonons, Joule heating, cooling to the substrate and thermal diffusion within the TES. 

P = PqP + ^Joulc + -Psubstratc iffusion 

(57) 

Most of these processes are relatively simple to model and discussed only briefly. Diffusion how- 
ever is more involved and needs to be considered if it is necessary to properly describe normal- 
superconducting phase separation in the TES. 



1. Phonon Heat 



The process of phonon - quasiparticle down conversion was described previously in Section [TV] as 
a process in which phonons energy is removed from the crystal. This process results in a population 
of quasiparticles (broken Cooper pairs) which diffuses through the aluminum films. In the CDMS 
detectors, these aluminum films are coupled to the tungsten TESs. Due to the large number of 
quasiparticles (the aluminum quasiparticle gap is ^100 /ieV, small compared to the 10-100 keV 
gamma energy) the details of where individual phonons interact in the film average out and we can 
use average quantities. In the CDMS detectors, the measured efficiency of quasiparticles reaching 
the TES is eqp ~10%. The remaining ^90% are trapped in the aluminum (details are discussed in 
references [S^JIICT] . 

^QP = eqp-Pph onon 

(58) 



2. Joule Heating 

After determining the voltages at each node. Joule heating is calculated by 



Pj = ^. (59) 




FIG. 16: TES resistor interconnects as modeled using a finite element approximation. 



34 



3. Substrate Cooling 

Power into the substrate is given by 

Ps = -S Vol (T^ES - T^Jubstratc) (60) 

where E is the electron-phonon coupUng parameter, Vol is the TES volume and n is an exponent 
that describes electron-phonon coupling. 

4. Diffusion Withm the TES 

Thermal diffusion is simplified by assuming that the QETs are thermally decoupled from each 
other due to poor conduction through the substrate. The process is then described by a 1-d diffusion 
process 

57dffi^ = DAT, (61) 

where 

AT = V • VT — ^ (62) 

and D is given by K/Cy where K is the thermal conductivity and Cy is the heat capacity per unit 
volume pi] . 

For non-boundary nodes, AT can be approximated by considering the midpoint and two adjoining 
points. To see this let us first consider just one dimension, specifically the term First we 

compute In the discrete approximation we can consider some point Ti and either its right 
neighbor or left neighbor T^+i or T^-i, where i represents the position of some node in the x 
direction. We obtain either 



or 



dT 
dx 

dT 
dx 



^ 5x ^^^^ 

right "-^ 



^ ' (64) 

loft 



where denotes the move from continuous to discrete space. Next we compute 

O Ol ^ dx I right dx I loft 

dx dx Sx 

{5xY ■ ^^^> 

Boundary nodes are computed differently. To conserve energy we impose the Neumann boundary 
condition 

VT|bo,„darv = 0- (66) 
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Making the discrete approximation is simple, we compute 



T, - T,_i = 



(67) 



for boundary nodes. 

T will be the temperature matrix of size n x m where n and m represents the nodes in the y and 
X directions respectively. It would be nice to convert the following procedures for computing AT 
into a matrix operation which considering Equations |65| and |67| is fortunately is quite simple. The 
matrix A^; is of size m x m and allows discrete calculations of AT. Specifically, 



(68) 



where A,, has the form 



/-I 
1 




1 

-2 1 
1 -2 1 









\ 






••• 1 -2 1 
\ 0---0 1 -1 / 



(69) 



This form is equivalent to Equations 65 and 67 and that the signs are correct for the boundary 
nodes to ensure heat moves from hot to cold nodes. 
The diffusion algorithm is then given by 



T„ +D 6t 



TA, 
{Sxf 



The algorithm should be stable as long as the Courant-Friedrichs-Lewy condition 

(Sxr 



St < 



2D 



(70) 



(71) 



is satisfied [521 [S3] but the Joule heating and substrate cooling terms result in a less stable algorithm 
and the stability is observed to become poorer as the TES cell size is reduced. The onset of instability 
is due to local temperature deviations due to small numerical imprecision. These deviations can be 
recognized by comparing the maximum absolute temperature difference in adjacent cells T^iff adj = 
max(|(ri-|_i — Til) to the maximal temperature difference in the TES T^iff^max = max(T) — min(r). 
The number of time steps is then scaled up by Tdiff^adj /2diff,max x N/ f, where N is the number of 
cells along the TES and / is a factor ^ 10 to reduce the instability caused by the inclusion of Joule 
heating and substrate cooling. 



D. Parsing Phonons Into The TESs 



One issue that needs to be considered before running a TES simulator is the parsing of phonons 
into the multiple QETs which are wired in parallel. It is generally more computationally efficient to 
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ignore the exact location with and without aluminum coverage and to use an average approximation. 
In this case, we need only determine the physical region of each QET structure. This can be 
approximated by assuming each QET to have the same area. For a TES channel with radial and 
angular extent R2, Qi and 02 and containing N discrete circuit elements, the area of each 
element is simply 

« = ^ (^2 - Rl) (02 - 61) . (72) 
For the first ring of TES elements we could also consider the area of each element to be 

1 sn 

a = -59[{Ri + drf - Rj] ^ —(2R^Sr + (^r^), (73) 

where and SB are the elements' size in the radial and angular directions. Equating the two 
descriptions of area we get 

^(i?^ - i??)(e2 - 61) = m [r^St + \br^\^ . (74) 

We can now impose a proportionality between the element's linear dimensions in the radial and 
theta direction (measured at r — Ri + Sr/2) then we get the relationship 5r = c(i?i + 6r/2)S9, 
where c is a proportionality constant. This can be substituted into Equation |74| to eliminate SO and 
we obtain 

N Sr"" + 2RiN 5r' - y (i?2 - RI){^2 - 61) 6r - y (i?2 " ■R?)(02 - ei)i?i = 0. (75) 

This equation can be solved for 6r and provides one real root. 

If (5r > i?2 — -Ri then the lower value R2 — Ri is of course selected otherwise a 5r for each ring is 
proposed to be equal and rounded off to (i?2 — i?i)/round((i?2 — Ri)/5r). The round is included 
to ensure an integer number of rings. The number of 9 divisions n for this ring is simply 

n = round((i?i + 5rf - RI{Q2 - ei)/2/a), (76) 

where the number of elements is further constrained by n > 1. 

This solution isn't really self consistent since we have imposed too many completing desires, 
namely that the radial and angular dimensions satisfy a proportionality condition and that there 
are an integer number of elements in a ring. The former can be relaxed, but the integer number of 
elements is mandatory. We can obtain a solution with dimensions close to the desired proportion- 
ality condition if we set 59 = (02 — Qi)/n and Sr = yj2a/59 + R\ — Ri. After this first ring is built 
up the procedure repeats, with i?i ^ i?i — 5r until the entire channel area has been segmented. 

E. Numerical Constants for TES Simulation 



Table III lists several constants related to the TES thermodynamic properties and quasiparticle 
collection efficiency from aluminum to tungsten. 
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TABLE III: Physical constants for tungsten TES and aluminum fin simulation, from reference [57] ■ 



name 


symbol 


value 


units 


n 


electron-phonon coupling exponent 


5 


n/a 


S 


electron-phonon coupling constant 


4.8x10** 


W m-3 K-'^ 


c 


heat capacity 


37 


J 


D 


diffusion constant 


4x10"" 


rr? s 


eqp 


QP detection efficiency 


10% 


n/a 



VIII. FINAL REMARKS 

This paper has covered many physics and Monte Carlo topics in cryogenic radiation-detectors 
that utilize phonon and ionization readout. The bulk response is different for gamma and neutron 
interactions, the former producing a larger ratio or ionization energy. Propagation of both the 
phonons and charge carriers is anisotropic, the phonons transport dependent on dispersion relations 
and the charge carrier transport dependent on mass tensors. Phonon interactions are complicated at 
the surface by electron-phonon downconversion requiring additional Monte Carlo effort. The TES 
readout simulation is relatively straightforward relying on basic thermal and electrical processes. 
Various numerical techniques including PDF sampling and tricks for improving efficiency have been 
discussed. 

Monte Carlo of detectors is a valuable tool for detector design, characterization and data analysis. 
Often counter intuitive results can be quickly discovered in Monte Carlo which aids in prioritizing 
laboratory R&D efforts. Furthermore, Monte Carlo can help to explain features in data, or at least 
rule out models. 

Detailed Monte Carlo studies in CDMS [li [371 US IM'-'SS; and EDELWEISS [13111] detectors 
can be found in the literature. These studies include discussion of quantities which are relevant in 
the tuning and their effects on phonon pulse shape, charge transport, TES phase-separation and 
TES noise. 

With the power of Monte Carlo and increased use of low temperature cryogenic detectors, a 
general package is becoming more appealing. The GEANT4 collaboration is beginning to implement 
phonon and charge transport modules into the GEANT4 toolkit and will expand the benefits of 
this research to others outside the cryogenic detector community |69| . 
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